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Abstract 

In these lectures we describe the use of Monte Carlo simulations in understanding the role of 
tunneling events, instantons, in a quantum mechanical toy model. We study, in particular, a 
variety of methods that have been used in the QCD context, such as Monte Carlo simulations of 
the partition function, cooling and heating, the random and interacting instanton liquid model, 
and numerical simulations of non-Gaussian corrections to the semi-classical approximation. 
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I. INTRODUCTION 



We consider a non-relativistic particle moving in a potential V(x). The Hamiltonian of 
the system is given by 

We can rescale x and t such that 2m = A = 1. We will also use h = 1. This means that the 
system is characterized by just one dimensionless parameter, rj. The potential V(x) with 
i] = 1.4 is shown in Fig. [1^. The physics of this system is easy to understand. Classically, 
there are two degenerate minima at x — ±77. Quantum mechanically, the two states can 
mix. If the potential barrier is very high, rj — > 00, the wave functions of the ground state 
and the first excited state are approximately given by 

Mx) = -j=(il>-(x)±il> + (x)), (2) 

where ip±(x) are the ground state wave functions in the left and right minimum of the 
potential. The energy splitting between the ground state and the first excited state is 
exponentially small. The WKB approximation gives 
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AE = E 1 -E = d—±wexp(-So), (3) 



where u = Arj and So = m 2 u 3 /(12X) = 4r) 3 /3. 

Applications of the WKB method and of instantons to t 
discussed in many reviews and text books □□□□□□□ 
present another review on the subject in these lecture notes. Instead, we will use the double 
well potential to illustrate a number of numerical methods that have proven useful in the 
context of QCD and other gauge theories. 



oe double well potential are 
. It is not our intention to 



II. EXACT DIAGONALIZATION 

The quantum mechanical problem defined by the Hamiltonian equ. (0) can be solved by 
determining the eigenvalues and eigenvectors of H. This can be achieved by choosing a basis 
and numerically diagonalizing the Hamilton operator in that basis. We have chosen a simple 
harmonic oscillator basis defined by the eigenstates of 

H = ^ + \mulx\ (4) 
2 



FIG. 1: a) Double well potential V(x) = (x 2 — rj 2 ) 2 for r\ = 1.4. We have indicated the position 
of the ground state and the first three excited state, b) Spectrum of the double well potential as a 
function of the parameter rj. In this figure we show the position of the first six states. We clearly 
observe that positive and negative parity states become degenerate as rj — > oo. 

The value of uq is arbitrary, but the truncation error of the eigenvalues computed in a finite 
basis will depend on the choice of ujq. In practice, however, this dependence is quite weak. 
The eigenstates of H satisfy H \n) = \n)uo(n + 1/2). The Hamiltonian of the anharmonic 
oscillator has a very simple structure in this basis. The only non-zero matrix elements are 

(n\H\n) = 3Ac 4 \(n + l) 2 + n 2 } + Bc 2 (2n + 1) + cu (n + 1/2) + c, (5) 



(n\H\n + 2} = Ac\±n + 6W(n + l)(n + 2) + Bc 2 J(n + l)(n + 2), (6) 



(n\H\n + 4) = C y(n + l)(n + 2)(n + 3)(n + 4), (7) 

as well as the corresponding hermitian conjugates. We have also denned A = 1, B = 
—2rj 2 — lJq/A, C = rf and c = 1/ '^/ujq. Note that both H and H conserve parity. We can 
decompose the matrix H nm = (n\H\m) into even and odd components, H = H even + H dd, 
such that the eigenvectors of H even and H odd have positive and negative parity, respectively. 

With the choice c^o = 4^ even modest basis sizes such as iV = 40 are sufficient in order 
to determine the first few eigenvectors very accurately. In Fig. \\)p we show the first six 
eigenvalues as a function of the parameter 7]. We clearly observe that as rj increases pairs 
of eigenvalues corresponding to even and odd eigenfunctions become almost degenerate. In 
this limit the eigenfunctions are of the form given in equ. (J2J). 
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III. QUANTUM MECHANICS ON A EUCLIDEAN LATTICE 



An alternative to the Hamiltonian formulation of the problem is the Feynman path inte- 
gral [9]. The path integral for the anharmonic oscillator is given by 

(xi \e~ iHtl \x ) = [ xih) ~ Xl Vx e *s s = f h dt ( \jA _ r x 2 _ 2^ _ (g) 
Jx(o)=x Jo V4 / 

In the following we shall consider the euclidean partition function 

Z(T) = J Vxe~ s \ S E = J'dr Qx 4 + (x 2 - ?? 2 ) 2 ) , (9) 

where (3 = 1/T is the inverse temperature. The partition function can be expressed in 
terms of the eigenvalues of the Hamiltonian, Z(T) = J2n ex P{~E n /T). In the following we 
shall study numerical simulations using a discretized euclidean action. For this purpose we 
discretize the euclidean time coordinate — ia, % — 1, . . . n. The discretized action is given 
by 

S = J2 {^(*< - *<-0 2 + ~ V 2 ) 2 } , (10) 

where Xi = x(ri). We shall consider periodic boundary conditions xq = x n . The discretized 
euclidean path integral is formally equivalent to the partition function of a statistical system 
of "spins" Xi arranged on a one-dimensional lattice. This statistical system can be studied 
using standard Monte-Carlo sampling methods. In the following we will simply use the 
Metropolis algorithm 11, 12]. The Metropolis method generates an ensemble of config- 



urations {xi}^ where i = 1, . . . , n labels the lattice points and k = 1, . . . , N con f labels the 
configurations. Quantum mechanical averages are computed by averaging observables over 
many configurations, 

N r 

(°) = N lim w~ £ ° {k) ( n ) 

N conf -*oc l\ con f k=1 

where is the value of the classical observable O in the configuration {xi}( k \ The 
configurations are generated using Metropolis updates {xj}^ — > {xj}^ +1 ^. The update 
consists of a sweep through the lattice during which a trial update x[ k+1 ^ = x^ + 8x is 
performed for every lattice site. Here, 5x is a random number. The trial update is accepted 
with probability 

P (xf ) -> xf = min {exp(-AS), 1} , (12) 

where AS is the change in the action equ. (JTUj). This ensures that the configurations {xj}^ 
are distributed according the "Boltzmann" distribution exp(— S). The distribution of Sx is 



arbitrary as long as the trial update is micro-reversible, i. e. is equally likely to change x[ 
to X} 1 ' and back. The initial configuration is also arbitrary. In order to study equilibration 
it is often useful to compare an ordered (cold) start with {xi}^ = {r/} to a disordered (hot) 
start {xi}^ = {rj}, where r, is a random variable. 

The advantage of the Metropolis algorithm is its simplicity and robustness. The only 
parameter to adjust is the distribution of 5x. We typically take 5x to be a Gaussian random 
number with the width of the distribution adjusted such that the average acceptance rate 
for the trial updates is around 50%. Fluctuations of O provide an estimate in the error of 
(O). We have 

A(0) = , {U L W . (13) 

\ -''con/ 

This requires some care, because the error estimate is based on the assumption that the 
configurations are statistically independent. In practice this can be monitored by computing 
the auto-correlation "time" in successive measurements 0({xi}^). The auto-correlation 
time of different observables can be very different. For example, successive measurements of 
the total energy decorrelate very quickly, but measurements of the topological charge have 
a much longer correlation time. 

The energy eigenvalues and wave functions of the quantum mechanical problem can be 
obtained from the euclidean correlation functions 

n(r) = (O(0)O(t)>. (14) 

Here, 0(t) is an operator that can be constructed from the variables x(t), e.g. some 
integer power 0(r) = x(r) n . The euclidean correlation functions are related to the quantum 
mechanical states via spectral representations. The spectral representation is obtained by 
inserting a complete set of states into the expectation value equ. (JHJ). We find 

n(r) = £ |(0|O(0)|n)| 2 exp(-(K - E )r), (15) 

n 

where E n is the energy of the state \n) and |0) is the ground state of the system. We can 
write this as 

II(t) = J dE p(E) exp(-(E - E q )t). (16) 

In the case we are studying here there are only bound states and the spectral function p(E) 
is a sum of delta-functions. Equ. (|T3jl shows that the euclidean correlation function is easy 
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FIG. 2: Typical euclidean path obtained in a Monte Carlo simulation of the discretized euclidean 
action of the double well potential for rj = 1.4. The lattice spacing in the euclidean time direction 
is a = 0.05 and the total number of lattice points is N T = 800. The green curve shows the 
corresponding smooth path obtained by running 100 cooling sweeps on the original path. 

to construct once the energy eigenvalues and eigenfunctions are known. This fact was used 
in order to calculate the solid lines shown in Figs. (J3J). The inverse problem is well denned in 
principle, but numerically much harder. In the following we will concentrate on extracting 
just the first few energy levels. A technique that can be used on order determine the spectral 
function from euclidean correlation functions is the maximum entropy image reconstruction 



method, see 



|l3, 



The Monte Carlo method is very useful in calculating expectation values in quantum or 
statistical mechanics. However, the Monte Carlo method does not directly give the partition 
function or the free energy. In principle one can reconstruct the free energy from the energy 
eigenvalues but this is not very practical since, as we just mentioned, it is hard to compute 
the full spectrum. A very effective method for computing the free energy is the adiabatic 
switching technique. The idea is to start from a reference system for which the free energy 
is known and calculate the free energy difference to the real system using Monte Carlo 
methods. 

For this purpose we write the action as S a = S + aAS where S is the full action, So 
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is the action of the reference system, AS" is defined by AS = S — So, and a is a coupling 
constant. The action S a interpolates between the real and the reference system. Integrating 
the relation d log Z (a) /(da) = —(AS) a we find 

log(Z(a = 1)) = log(Z(a = 0)) - f da 1 (AS)* , (17) 

where (.) a is an expectation value calculated using the action S a . In the present case it is 
natural to use the harmonic oscillator as a reference system. In that case 

z(a = 0) = £gcp( _^ )= exp(-/3W2) 
V 1 - exp(-(3uj ) 

where c^o is the oscillator constant. Note that the free energy of the anharmonic oscillator 
should be independent of u)q. The integral over the coupling constant a can easily be calcu- 
lated in Monte Carlo simulations by slowly changing a from to 1 during the simulation. 
In order to estimate systematic errors due to incomplete equilibration it is useful to repeat 
the calculation with a changing from 1 to and study possible hysteresis effects. 



IV. NUMERICAL RESULTS 



Numerical results from Monte Carlo simulations of the euclidean path integral are shown 
in Figs.EHU The numerical data were obtained using the program qm. f or which is described 
in more detail in the appendix. A typical path that appears in the Monte Carlo simulation 
is shown in Fig. El The figure clearly shows that there are two characteristic time scales in 
the problem. On short time scales the motion is controlled by the oscillation time t osc ~ 
u^ 1 ~ (4r/) _1 . For large r the system is governed by the tunneling time r tun ~ exp(— 4rj 3 /3). 
In order to perform reliable simulations we have to make sure that the lattice spacing a is 
small compared to t osc and that the total length of the lattice Na is much larger than the 
tunneling time 

a < t osc , r tun < Na. (19) 

A typical choice of parameters for the case r] = 1.4 is a number of lattice points N = 800, a 
lattice spacing a = 0.05 and a number of Metropolis sweeps N con f = 10 5 . 

Fig-Olshows the distribution of Xi obtained in the Monte Carlo simulation compared to the 
square of the ground state wave function computed by the diagonalization method discussed 
in section |H] As rj is increased and the potential barrier becomes larger the tunneling time 
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FIG. 3: Probability distribution |-0(x)| 2 in the double well potential for r/ = 1.4. The solid line 
shows the "exact" numerical result obtained by diagonalizing the Hamiltonian in an oscillator basis 
whereas the histogram shows the distribution of x for an ensemble of euclidean paths. 

increases exponentially and the number of configurations needed to reproduce the correct 

wave function becomes very large. 

Fig. 0] shows the correlation functions of the operators x, x 2 and x 3 . The solid lines show 

the result obtained using the spectral representation equ. (j!5|) together with the eigenvalues 

and eigenfunctions determined by numerical diagonalization of the Hamiltonian. The data 

points show the results from the Monte Carlo simulation. There is a small systematic 

disagreement for small r which is related to discretization errors but the overall agreement is 

excellent. Energy levels and matrix elements can be obtained from the logarithmic derivative 

of the correlation function, 

rflogn(r) = gnOgn ~ Eo)|(0|Q(0)|n)| 2 exp(-(E w - E )r) 
{T> " dr " Zn\(0\O(0)\n)\ 2 eM-(E n -E )T) ' 1 > 

In the limit r — > oo the function C (r) converges to the energy splitting between the ground 

state and the first excited state that has a non- vanishing transition amplitude (0|O(0)|n). 

Because of parity invariance, O = x n connects the groundstate to parity even/odd levels for 

n even/odd. Since the first excited state is parity odd we have 

lim — log(x(r)x(0)> = lim — log(x 3 (r)x 3 (0)) = E 1 - E . (21) 




FIG. 4: Fig. a shows the correlation functions (O(O)O(r)) in the double well potential for rj = 
1.4 and O = x,x 2 ,x 3 . The solid lines are "exact" numerical results obtained by diagonalizing 
the Hamiltonian in an oscillator basis whereas the data point were obtained from Monte Carlo 
simulations with a = 0.05 and N T = 800. Fig. b shows the logarithmic derivative of the correlators 
in Fig. a. In the case of the (x 2 (0)x 2 (r)) we subtracted the constant contribution. 

For even powers of x the situation is more complicated because the correlator has a constant 
term |(0|x 2n |0)| 2 . After subtracting the constant part, the logarithmic derivative of the 
correlation function of even powers of x tends to (E 2 — E ). Numerical results are shown 
in Fig. We observe that the logarithmic derivative of (x(r)x(O)) converges very rapidly 
to AEi = (Ei — Eq). The numerical results for AE 2 = (E 2 — E Q ) have large uncertainties. 
These uncertainties are related to the fact that the correlator (a; 2 (r)a;(0)) is dominated by 
the subtraction constant (x 2 ) 2 . The logarithmic derivative of the (x 3 (r)x 3 (0)) correlator also 
tends to AEi, but receives larger contributions from excited states. This feature can be used 
in order to extract the energies of higher states. The idea is very simple. From the matrix 
elements C\ = (0|x|l) and d\ = (0|a; 3 |l) we can determine a new operator O — x/ci — x 3 /di 
that does not couple to the first excited state. This operator predominantly couples to 
the third excited state. Repeating this procedure we can determine the energies of higher 
excited states. The problem is that correlation functions of higher powers of x are more and 
more noisy. As a result, finding the energies of highly excited states is very hard, even in 
the simple problem considered here. 
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FIG. 5: Free energy F = —T\og(Z) of the anharmonic oscillator as a function of the temperature 
T = 1/(5 with (5 = na. The solid line was calculated using the spectrum of the Hamiltonian. The 
data points were obtained using Monte Carlo calculations and the adiabatic switching method. 

In Fig. El we show Monte Carlo results for the partition function compared to "exact" 
results based on the spectrum of the anharmonic oscillator obtained in Sect. |Hj The Monte 
Carlo results agree with the direct calculations but the Monte Carlo method is effectively 
limited to a small range of temperatures. If the temperature is very small the partition 
function is dominated by the ground state contribution. In that case, it is much more 
efficient to compute the ground state energy directly by measuring the expectation value 
of the Hamiltonian, Eq = (H) , with H = x 2 /A + V(x). There is one subtlety with this 
approach: If a naive one-sided discretization of the time derivative is used then the continuum 
limit of the expectation value of the kinetic energy diverges. This problem can be addressed 
by using an improved discretization of the kinetic energy or by using the Virial theorem. 
The Virial theorem implies that 

{H) = {T + V) = { X -V' + V). (22) 

At high temperature more and more states contribute. The main difficulty with the Monte 
Carlo approach in this regime is that discretization errors have to be carefully monitored. 
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FIG. 6: Same as Fig. 0] but the correlation functions are evaluated from cooled Monte Carlo 
configurations. The number of cooling sweeps is N coo i = 200. 

V. EXTRACTING THE INSTANTON CONTENT USING COOLING 



From Fig. |21 we can clearly see that for this particular choice of the parameter r\ a typical 
path contains two components, one related to quantum fluctuations with frequency u, and 
one related to tunneling events, instantons. In the continuum limit the instanton solution 
can be found from the classical equation of motion 

S 



Sx(t 



-S i 







mi = V'(x). 



(23) 



The solution which satisfies the boundary condition x(r — > ±oo) = ±77 is given by 



XjlT) 



?7tanh 



(24) 



where u = 4i] and tq is the "location" of the instanton. The anti-instanton solution is simply 
given by Xa(t) = —xi(t). The classical action of the instanton is 



(25) 



The tunneling rate n I+ A = N i+ a/P is exponentially small, n I+ A ~ exp(— Sq). In order to 
determine the pre-factor one has to study small fluctuations around the instanton solution. 
This calculation has been carried out to next-to-leading order in the semi-classical expansion. 
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The result 

^ = 8,V^ exp (_ So _ZIi_). (26) 

The tunneling events can be studied in more detail after removing short distance fluctuations. 
A well known method for doing this is "cooling" 17]. In the cooling method we only 



accept Metropolis updates that lower the action. This will drive the system towards the 
nearest classical solution. Since instantons are classical solutions, cooling can be used to 
study the instanton content of a quantum configurations. This is clearly seen in Fig. EJ The 
black line is the original, quantum, configuration. The green line is the same configuration 
after 200 cooling sweeps. It is easy to check that this configuration is very close to a linear 
superposition of independent tunneling events. For this purpose we can extract the instanton 
and anti-instanton locations from the zero crossings and compare the cooled configuration 
to the simple "sum ansatz" 



X .s 



( r ) = v \ Qi tann 



(27) 



where Qi = ±1 is the topological charge of the instanton. The most important question is to 
what extent physical observables in the cooled configurations resemble those in the original 
configurations. This provides a measure of the importance of instantons in the double well 
potential. In Fig. El we show correlation functions measured in the cooled configurations. 
These results should be compared with the full correlation functions shown in Fig. HJ We 
observe that the correlation functions are quite different. Short distance fluctuations elim- 
inated by cooling obviously play an important role. We observe, however, that the level 
splitting between the ground state and the first excited state is clearly dominated by semi- 
classical configurations. The logarithmic derivative of both the (x(0)x(t)) and (x 3 (0)x 3 (r)) 
correlation functions is very well reproduced in the cooled configurations. 



VI. THE DENSITY OF INSTANTONS 



The cooling method can also be used in order to get an estimate of the total density of 
instantons and anti-instantons. While the net topological charge, the number of instantons 
minus the number of anti-instantons, is unambiguously defined the same is not true for the 
total number of topological objects. There is no clear distinction between a large quantum 
fluctuation and a very close instanton-anti-instanton pair. In the cooling method this is 

12 




FIG. 7: Instanton density and instanton action as a function of the number of cooling sweeps for 
different values of the parameter ry. The solid and dashed green lines in Fig. a shows the semi- 
classical instanton density at one and two- loop order. The solid line in Fig. b shows the classical 
instanton action. 

reflected by the fact that the number of instantons, extracted from the number of zero 
crossings in the cooled configuration, depends on the number of cooling sweeps. 

It is clear, however, that the instanton density should be well denned in the semi- 
classical limit. In this limit there is an exponentially large separation of scales between 
the tunneling time r tun and the scale of ordinary quantum fluctuations r osc . This separa- 
tion of scales can also be exploited in the cooling method. Cooling is a local algorithm 
which implies that it takes on the order of r/a cooling sweeps in order to affect coher- 
ent structures that exist at a scale r. We expect that the number of instantons mea- 
sured using the cooling method is approximately given by the sum of two exponentials, 
Niiricooi) = N osc exp(-n cool a/r osc ) + N tun exp(-n cool a/r tun ). The first exponential describes 
the disappearance of quantum fluctuations on a time scale t osc and the second exponential 
reflects instanton- anti- instanton annihilation occurring on a time scale T tun . 

Numerical results for Ni(n coo i) are shown in Fig. [3 We observe that the data is consistent 
with the presence of two distinct time scales and that the description in terms of two expo- 
nentials becomes better as the semi-classical limit rj — > oo is approached. We also note that 
after the quantum noise has disappeared the instanton density is close to the semi-classical 
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FIG. 8: Instanton density as a function of the parameter r]. The blue symbols show the instanton 
density extracted from Monte Carlo configurations after 10 cooling sweeps. The red symbols show 
the results of a Monte Carlo calculation of non-Gaussian effects. The green lines show the semi- 
classical instanton density at one and two-loop order. The black line shows AE/2 where AE is 
the splitting between the ground state and the first excited state. 

result equ. (|26j) . A more detailed comparison is shown in Fig. [SJ In this figure we show the 
instanton density after 10 cooling sweeps, the one and two-loop semi-classical result, as well 
as the level spacing between the ground state and the first excited state. 

We observe that for r\ > 1.2, corresponding to a classical instanton action So > 2, the 
number of instantons extracted using the cooling method agrees very well with the semi- 
classical approximation. We also note that the two-loop result is a clear improvement over 
the one-loop approximation for classical actions as small as Sq ~ 1. Finally, we observe that 
the instanton density is close to the level splitting even in the regime where So is less than 
one. 

In Fig. [7| we also show a Monte Carlo calculation of the instanton density on a small 
lattice. The idea is very simple. The one-loop calculation of the tunneling rate is based on 
expanding the action around the classical path to quadratic order 

5 2 S 



S = S + - J dr 5x(t) 



5x 2 



5x(t) + . . . , (28) 

xj(t) 



14 



-3 







2 



3 



4 



5 



x 



FIG. 9: Quantum mechanical paths which appear in a Monte-Carlo calculation of the one-instanton 
partition function in the double well potential. The calculation involves adiabatic switching be- 
tween the Gaussian effective potential and the full potential. The smooth curves are the initial 
configurations in the zero and one-instanton sector. The Monte Carlo updates in the one-instanton 
sector involve a constraint which keeps the instanton location fixed. 

where Sx(t) = x(t)—xj(t). As in Sect. lIIII we can introduce a new action S a that interpolates 
between the full action and the Gaussian approximation, S a = S gauss + aAS with AS = 
S — S gauss . The exact quantum weight of an instanton can be determined by integrating 
over the coupling constant a. We have 



where (.)^ is an expectation value in the n-instanton sector at coupling a. The method 
is illustrated in Fig. |H1 The figure shows typical paths that contribute to (AS) in the zero 
and one-instanton sector. The resulting estimate of the instanton density is also shown in 
Fig. |H| The Monte Carlo results show that the instanton density is reduced compared to 
the one-loop estimate. For classical instanton actions So > 3 the result is in agreement 
with the two-loop estimate and the cooling calculation. It is hard to push the Monte Carlo 
calculation to instanton actions Sq < 3 because transitions between the zero and two (four, 
six, . . .) sector become too frequent. 




(29) 
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FIG. 10: Same as Fig. |1] but the correlation functions are evaluated from a random instanton 
configuration. 

VII. THE INSTANTON LIQUID MODEL 

Given the success of the semi-classical approximation in predicting the splitting between 
the ground state and the first excited state it seems natural to study the correlation functions 
in the semi-classical approximation in more detail. We begin by considering the contribution 
from the classical path only. In this case the partition function is given by 

^£^(n/4« (30) 

Here, ti^ua are the number of instanton and anti-instantons, Tj are the (anti) instanton 
positions, and So is the classical action. In the next section we will discuss the problem of 
choosing the correct path for a multi-instanton configuration in more detail. The simplest 
choice is the sum ansatz given in equ. (J27j) . The coordinate correlation function is given by 

n rf (r) = (x d {0)x cl {r)), (31) 

where (.) denotes an ensemble average over the collective coordinates r,. The distribution 
of collective coordinates is controlled by the partition function equ. (J3(Jj) . The simplest 
approximation is to ignore the interaction between instantons. In this case the action is 
S = (nj + riA)S and the distribution of collective coordinates is random. This is known as 
the instanton gas model or the random instanton approximation. 
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FIG. 11: Gaussian effective potential for small fluctuations around a single instanton path centered 
at r = 0. 

Correlation functions in a random instanton gas are shown in Fig. El We note that 
the result is very similar to the cooled correlation functions shown in Fig. H This is in 
agreement with our earlier observation that the cooled configurations are very close to a 
simple superposition of instantons. Like the the cooling calculation the random instanton 
gas reproduces the splitting between the ground state and the first excited state, but it does 
not give a good description of other aspects of the correlation functions. 

It is clear that the main feature that is missing from the ensemble of classical paths 
is quantum fluctuations. Quantum fluctuations appear at next order in the semi-classical 
approximation. We already noted that quantum fluctuations determine the pre-exponential 
factor in the tunneling rate, see equ. (J2BJ). We can write the path as x(t) = x c /(r) + 8x(t) 
where x c i(r) is the classical path and 8x(t) is the fluctuating part. To second order in 
5x the action is given by equ. (|28|). For a single instanton it is possible to determine the 
propagator (Sx(0)Sx(t)) analytically, see equ. (39) in Q|. For an ensemble of instantons 
we can approximate the propagator as a sum of contributions due to individual instantons. 



2l|. 



This is the procedure that is used in the QCD calculations described in 

Alternatively, we can determine the correlation function numerically, using the "heating" 
method. As the name suggests, this is essentially the inverse of the cooling method. We begin 
from a classical path and determine the Gaussian effective potential for small fluctuations 
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FIG. 12: Typical random instanton configuration and the same configuration with Gaussian 
fluctuations. The noisy path was generated using 10 heating sweeps in the Gaussian potential 
around the classical path. This figure should be compared with Fig. 

around the path. For a single instanton, the action is given by 
S = Jdr ^5x 2 (r) + V 

see Fig. This action has one zero mode 5x(r) = —S 1/z dx c i(T — Ti)/(drj) which cor- 
responds to translations of the instanton solution. We can eliminate the corresponding 
non-Gaussian fluctuations by imposing a constraint on the location of the instanton. Using 
the simple identity 

1 = JdnSixin)) \x(rj)\ (33) 

we see that the corresponding Jacobian is the velocity x(ti). We can now perform Monte 
Carlo calculations using the Gaussian action for a multi-instanton configuration. The 
method is illustrated in Fig. The black line shows the classical path and the green 
path is the same path with Gaussian fluctuations included. Clearly, this path looks very 
similar to the full quantum path shown in Fig. |21 There are still some differences, how- 
ever. We notice, in particular, that the fluctuations around the minima of the potential are 
not completely symmetric. This is related to non-Gaussian effects. We also observe that 
the heated random instanton path lacks large excursions from the minima of the potential 
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2cosh 2 (2?7(T - n)) 



6x 2 (t) 



(32) 




FIG. 13: Same as Fig.0]but the correlation functions are evaluated in a random instanton ensemble 
with Gaussian fluctuations. 

that do not lead to a tunneling event. These effects are due to a combination of instanton 
interactions and large non-Gaussian effects. 

Correlation functions in the random instanton configurations with Gaussian fluctuations 
included are shown in Fig. ED We observe that the correlation functions are in much better 
agreement with the exact results than the correlators obtained from the classical path only. 
We also see that the correlators not only describe the splitting between the ground state 
and the first excited state but also provide a reasonable description of the second excited 
state. 



VIII. INSTANTON INTERACTIONS 

Another feature that is missing from the random instanton ensemble is the correlation 
between tunneling events due to the interaction between instantons. In QCD instanton 
interactions, in particular those mediated by fermions, are very important and lead to qual- 
itative changes in the instanton ensemble. In the quantum mechanical model studied here 
the interaction between instantons is short range and only leads to relatively small effects. 
These effects can nevertheless be clearly identified in very accurate calculations. We refer to 
[l] for a discussion of the contribution of instanton-anti-instanton pairs to the ground state 
energy. 
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FIG. 14: Instanton-anti-instanton interaction in units of So as a function of the instanton-anti- 
instanton separation. The solid line shows the result in the sum ansatz. The triangles show the 
same data plotted as a function of the zero crossing distance. The streamline interaction is shown 
as the circles and the squares show the effective interaction extracted from the cooled intanton- 
anti-instanton distribution shown in Fig. 1161 

The simplest method for studying the instanton-anti-instanton interaction is to construct 
a trial function and compute its action. For the sum ansatz given in equ. (|27jl the result 
is shown as the solid line in Fig. EI Asymptotically, the action is given by Sia(tia) — 
26*0(1 — 6 exp(— r]TjA) + • • •) where tja = |tj — ta\ is the instanton-anti-instanton separation. 
In the opposite limit, tja — * 0, the instanton and anti-instanton annihilate and the action 
tends to zero. It is clear, however, that in this limit the sum ansatz is at best an approximate 
solution to the classical equation of motion, and it is not obvious how the path should be 
chosen. 

The best way to deal with this problem is the "streamline" or "valley" method 22 , |3| . 
The method is based on the observation that in the space of all instanton-anti-instanton 
paths there is one almost flat direction along which the action slowly varies between 2Sq and 
0. All other directions correspond to perturbative fluctuations. We can force the instanton- 
anti-instanton path to descend along the almost flat direction by adding a constraint 

^ = £(A) / dr (x(r) - x A (r)) (34) 
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FIG. 15: Solution of the streamline equation for an instanton-anti-instanton pair. Figure a 
shows the streamline path and the Fig. b shows the action density. The paths correspond to 
S/S = 2.0,1.8, ... ,0.2,0.1. 

to the classical action. Here, A labels the different instanton-anti-instanton paths along 
the streamline and £(A) is a Lagrange multiplier. We find the streamline configuration by 
starting from a well separated IA pair and letting the system evolve using the method of 
steepest descent. This means that we have to solve 



with the boundary condition that x\=q(t) ~ x sum (r) corresponds to a well separated 
instanton-anti-instanton pair. Note that £(A) is an arbitrary function that reflects the 
reparametrization invariance of the streamline solution. A sequence of paths obtained 
by solving equ. (J3H|) numerically is shown in Fig. ED We also show the action density 
s = x 2 /4 + V(x). We can see clearly how the two localized solutions merge and eventually 
disappear as the configuration progresses down the valley. 

There is no unique way to parametrize the streamline path and extract the instanton- 
anti-instanton action as a function of the separation between the tunneling events. The 
simplest possibility is to used the distance between the zero crossings t z . This definition 
has the advantage of being very easy to use, but it prevents us from exploring the part of 
the streamline trajectory where the instanton and anti-instanton are so close that the path 
never crosses zero. In Fig. El we compare results for Sia{t z ) obtained from the sum ansatz 




(35) 
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FIG. 16: Distribution of instanton-anti-instanton separations after 10 cooling sweeps. The sepa- 
ration was extracted from the distance of the zero-crossings in the cooled configuration. 

and the streamline solution. We observe that for instanton separations t z > 0.3 the results 
are very similar. We also note, however, that for r z < 0.6 the zero crossing distance is quite 
different from the parameter n — ta that appears in the sum ansatz. 

One can show that the ambiguities that arise in trying to define the instanton-anti- 
instanton interaction at short distance correspond to similar ambiguities that arise in the 
perturbative expansion because and are related to the factorial growth of higher order terms 
in the expansion. Only the sum of the perturbative and the instanton contribution is well 
defined and leads to unique predictions for the groundstate energy. In these lectures we 
shall not discuss this problem any further. Instead, we will study the question whether the 
full quantum configurations contain evidence of the correlations between instantons that the 
classical instanton-anti-instanton interaction implies. 

In Fig. EI] we show a histogram of the instanton-anti-instanton separation determined 
in Monte Carlo simulations of the quantum mechanical partition function. The data were 
obtained by measuring the zero crossing distance after 10 cooling sweeps. For comparison 
we also show the I A distribution in the random instanton gas. We observe that there is an 
enhancement of close I A pairs which corresponds to an attractive instanton-anti-instanton 
interaction. The exact magnitude of this enhancement depends sensitively on the number of 
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FIG. 17: Typical instanton configuration in an interacting instanton calculation. The figure shows 
the location x of the first 10 instantons (blue) and anti-instantons (red) over a period of 3000 
configurations. 

cooling sweeps. As emphasized in the previous paragraph, this is not necessarily a shortcom- 
ing of the cooling method. We can try to translate the enhancement in the I A distribution 
into an effective interaction using the classical relation ti(tia) ~ ^o(tta) exp(— Sia(tia))- 
Here, n{r IA ) is the I A distribution, n (rjA) is the distribution in the random theory and 
Sia{ t ia) is the instanton-anti- instanton interaction. The result is also shown in Fig. EI We 
observe that the interaction extracted from the I A distribution is significantly weaker than 
the classical result. This may imply that the full quantum interaction is weaker than the 
classical result, or that too many close pairs are lost during cooling. This question can be 
studied in more detail using the methods discussed in Sect. IVII 

Finally, we address the question how to include correlations between tunneling events in 
the instanton calculation. For this purpose we include the instanton-anti-instanton inter- 
action in the instanton liquid partition function equ. (jHOJ)- In this context we again have 
to address the problem of close instanton-anti-instanton pairs. The simplest approach is to 
add a short range repulsive core which excludes configurations that are not semi-classical. 
The hard core interaction can be adjusted in order to reproduce the I A distribution found 
in the cooling calculation. In practice we have chosen S core (rj A ) = A c exp(— ti A /t c ) with 
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A c = 3 and r c = 0.3. In Fig. El we show a typical set of instanton and anti-instanton 
trajectories from an interacting instanton calculation. Correlations between instantons are 
clearly visible. In particular, we observe that a number of close instanton-anti-instanton 
pairs are formed. We have also studied correlation functions in the interacting instanton 
ensemble. We find that differences as compared to the random ensemble are rather subtle 
and a detailed study of non-Gaussian effects is necessary in order to establish the importance 
of instanton interactions. 

IX. SUMMARY 

In these lectures we presented Monte Carlo methods for studying the euclidean path 
integral in Quantum Mechanics. We also supply a set of computer codes written in fortran 
that were used to generate the data shown in the figures. We encourage the reader to play 
around with these programs in order to get a deeper appreciation of the path integral and 
of Monte Carlo methods. 

We should note that Monte Carlo calculations of the euclidean path integral are an 
extremely poor way to compute the spectrum or the correlation functions of the anharmonic 
oscillator. The code based on diagonalizing the Hamiltonian is both much faster and much 
more accurate than the Monte Carlo codes. The purpose of the Monte Carlo codes is 
entirely pedagogical. However, if we proceed from quantum mechanics to systems involving 
many more degrees of freedom, such as four- dimensional field theories, Hamiltonian methods 
become more and more impractical and Monte Carlo calculations based on the euclidean path 
integral provide the most efficient method for computing the spectrum and the correlation 
functions known to date. 

We also discussed Monte Carlo methods for studying the contribution of instantons to 
the euclidean path integral. In the case of the double well potential there is a parameter, 
77, which controls the instanton action Sq = Ar] 3 /3. If Sq ^> 1 then instantons are easily 
identified but the tunneling rate is small. If Sq ~ 1 then instantons are very abundant but 
it is hard to determine the instanton density precisely. We focused on the case Sq ~ 3 which 
is at the boundary of the semi-classical regime. Even though the expansion parameter 1/Sq 
is not very small the instanton density is still well determined and agrees with the level 
splitting. We also noticed, however, that non-Gaussian effects are important in this regime. 
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Ultimately, we are interested in the question to what extent these results can be gen- 
eralized to QCD. In QCD there is no free parameter that controls the applicability of the 
semi-classical expansion. Unlike the case of the double well potential, instantons in QCD 
can have any size. Asymptotic freedom implies that the action of small instantons is big, 
but the action of instantons with size p ~ ^qcd ls °f or der one. Nevertheless, lattice cal- 
culations support the idea that the tunnelin g de nsity is sizable, (N/V) ~ Aq CD ~ lfm~ 4 , 
and that instantons do not strongly overlap 24|. The typical instanton size is found to be 



p ~ (0.3 — 0.4) fm which implies a typical instanton action So ~ (5 — 10). The reason that 
the density is big even though the action is significantly larger than one is related to the 
fact that QCD instantons have many more collective coordinates, 12 (4 coordinates, 1 size, 7 
color angles) compared to just one in the case of the double well potential. As a consequence 
the pre-exponential factor in the tunneling rate is numerically large. 

There are some important differences between QCD and the double well potential. In a 
typical lattice QCD configuration quantum fluctuations of the gauge field are much bigger 
than the classical gauge fields associated with instantons. This implies that one cannot "see" 
instantons in the gauge configurations in the same way that one can immediately identify 
tunneling events in the quantum mechanical paths. Compare, for example, Fig. |21 with 
Fig. 1 in j^J. Only after some amount of cooling do instantons emerge from the quantum 
noise. On the other hand, fermions provide an important diagnostic tool in QCD that is not 
available in the simple bosonic model analyzed in these lectures. Instantons lead to localized 
chiral zero modes of the Dirac operators that can easily be identified even in noisy quantum 
configurations. 
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APPENDIX A: COMPUTER CODES 

All programs are written in standard fortran 77, have extensive comments and do not require 
any libraries. Some of the programs contain subroutine for generating random numbers or 
for diagonalizing matrices that were taken from Numerical Recipes (Numerical Recipes in 
Fortran, W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. D. Flannery, Cambridge 
University Press). 

1. qmdiag.for 

This programs computes the spectrum and the eigenfunctions of the anharmonic oscillator. 
The results are used in order to compute euclidean correlation functions. 
Input: fort. 05 

/ minimum of anharmonic oscillator potential V(x) = (x 2 — f 2 ) 2 

N dimension of basis used for diagonalizing H (choose N > 40) 

uo unperturbed oscillator frequency (choose ojq ~ 4/) 

Output: qmdiag.dat 



E n 


eigenvalue of Hamiltonian 




dipole matrix element (? n = (0 \x n) 2 


d n 


quadrupole matrix element d 2 n = \(0\x 2 \n)\ 2 




quadrupole matrix element d\ = (0 \x 3 \n) 2 


■ip(x) 


ground state wave function 


(x(O)x(r)} 


euclidean correlation function, also for x 2 and x 3 


d\ogU/(dr) 


log derivative of n(r) = (x(O)x(t)) 


m 


partition function 
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2. qm.for 

This programs computes correlation functions of the anharmonic oscillator using Monte 
Carlo simulations on a euclidean lattice. 
Input: fort. 05 

/ minimum of anharmonic oscillator potential V(x) = (x 2 — f 2 ) 2 

n number of lattice points in the euclidean time direction (n ~ 800) 

a lattice spacing (a ~ 0.05) 

ih ih = cold start xi = —f, ih = 1 hot start X{ = ran() 

n eq number of equilibration sweeps before first measurement (n eq ~ 100) 

n mc number of Monte Carlo sweeps {n mc ~ 10 5 ) 

5x width of Gaussian distribution used for Monte Carlo update x\ n ^ — > a;- n+1 - ) 

(6x ~ 0.5) 

n p number of points on which correlation functions are measured 

(xiX i+1 ), (xiX i+np ) (n p ~ 20) 
n mea number of measurements of the correlation function in a given Monte 

Carlo configuration Xi {n mea ~ 5) 
n pri number of Monte Carlo configurations between output of averages to 

output file (n pri ~ 100) 

Output: qm.dat 

S to t average total action per configuration 

V av ,T av average potential and kinetic energy 

(x n ) expectation value (x n ) (n = 1, . . . , 4) 

n(r) euclidean correlation function Il(r) = (O(O)O(r)) for O = x, x 2 , 

x 3 . Results are given in the format r, Il(r), All(r), d log(Il) /(dr), 
A[d log(n) / (dr)], where AIT(t) is the statistical error in n(r). 
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3. qms witch. for 

The program qmswicth.f or computes the free energy F = — Tlog(Z) of the anharmonic 
oscillator using the method of adiabatic switching between the harmonic and the anharmonic 
oscillator. The action is S a = Sq + a(S — So). The code switches from a = to a = 1 and 
then back to a — 0. Hysteresis effects are used in order to estimate errors from incomplete 
equilibration. Most input parameters are the same as in qm.for. Additional parameters are 
given below. 
Input: fort. 05 

uoq oscillator constant of the reference system (tu ~ 4/) 

n SW itch number of steps in adiabatic switching (n switc h ~ 20) 
Output: qmswitch.dat 

The output file contains many details of the adiabatic switching procedure. The final result 
for the free energy is given as F = F + 5F, where F is the free energy of the harmonic 
oscillator and 5F is the integral over a. We estimate the uncertainty in the final result 
as F ± AF(stat) ± AF(equ) ± AF(disc), where AF(stat) is the statistical error, AF(equ) 
is due to incomplete equilibration (hysteresis), and AF(disc) is due to discretizing the a 
integral. 

4. qmcool.for 

This programs is identical to qm.for except that expectation values are measured both in 
the original and in cooled configurations. We only specify additional input parameters. 
Input: fort. 05 

n s t number of Monte Carlo configurations between successive cooled config- 

urations. The number of cooled configurations is n con f/n st (n st ~ 20). 
n coo i number of cooling sweeps {n coo i ~ 50) 

Output: qmcool . dat 

Il(r) correlation functions are given in the same format as in qm.dat 

N[ + a total number of instantons extracted from number of zero crossings 

as a function of the number of cooling sweeps 
S t ot total action vs number of cooling sweeps 

S/N action per instanton. So is the continuum result for one instanton 
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5. qmidens.for 

The program qmidens.for calculates non-Gaussian corrections to the instanton density 
using adiabatic switching between the Gaussian action and the full action. The calculation 
is performed in both the zero and one-instanton-sector. The details of the adiabatic switching 
procedure are very similar to the method used in qmswitch.f or. Note that the total length 
of the euclidean time domain, (5 = na, cannot be chosen too large in order to suppress 
transitions between the one-instanton sector and the three, five, etc. instanton sector. Most 
input parameters are defined as in qm.for. 
Input: fort. 05 

^switch number of steps in adiabatic switching (n switc h ~ 20) 
Output: qmidens . dat 

The output file contains many details of the adiabatic switching procedure. The final result 
for the instanton density is compared to the Gaussian (one-loop) approximation. Note that 
the method breaks down if / is too small or (5 is too large. 
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6. rilm.for 

This program computes correlation functions of the anharmonic oscillator using a random 
ensemble of instantons. The multi-instanton configuration is constructed using the sum 
ansatz. Note that, in contrast to RILM calculations in QCD, the fields and correlation 
functions are computed on a lattice. 
Input: fort. 05 

/ minimum of anharmonic oscillator potential V(x) = (x 2 — f 2 ) 2 

n number of lattice points in the euclidean time direction (n ~ 800) 

a lattice spacing (a ~ 005) 

Nj + a number of instantons (has to be even). The program displays the one 



n 



n. 



n 



n 



mc 



rnea 



and two-loop result for the parameters (/, /3 = na). 
number of configurations (n mc ~ 10 3 ) 

number of points on which correlation functions are measured 

{xiX i+1 ), (xiX i+np ) (n p ~ 20) 

number of measurements of the correlation function in a given Monte 
Carlo configuration Xi (n mea ~ 5) 

number of Monte Carlo configurations between output of averages to 
output file {n pri ~ 100) 



Output: ri 

Stot 



lm. 



dat 



average total action per configuration 
average potential and kinetic energy, 
expectation value (x n ) (n = 1, . . . , 4) 



V T 



av 




H(r) 



euclidean correlation function Il(r) = (O(O)O(t)) for O = x, x 2 , 
x 3 . Results are given in the format r, n(r), An(r), d\og(Jl)/(dr), 
A[dlog(II)/(dr)]. 
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7. rilm_gauss.for 

This program generates the same random instanton ensemble as rilm.for but it also in- 
cludes Gaussian fluctuations around the classical path. This is done by performing a few 
heating sweeps in the Gaussian effective potential. Most input parameters are defined as in 
rilm.for. Additional input parameters are given below. 

Input: fort. 05 

n heat number of heating steps {n heat ~ 10) 

Sx coordinate update (Sx ~ 0.5) 

8. iilm.for 

This program computes correlation functions of the anharmonic oscillator using an inter- 
acting ensemble of instantons. The multi-instanton configuration is constructed using the 
sum ansatz. The configuration is discretized on a lattice and the total action is computed 
using the discretized lattice action. Very close instanton-anti-instanton pairs are excluded 
by adding an nearest neighbor interaction with a repulsive core. Most input parameters are 
the same as in rilm.for. Additional input parameters are 

Input: fort. 05 

r core range of hard core interaction (r core ~ 0.3) 

A core strength of hard core interaction (A core ~ 3.0) 

dz average position update (dz ~ 1) 



32 



